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We derive two different sets of rate equations for sympathetic cooling of 
harmonically trapped Bosons or Fermions. The rate equations are obtained 
from a master equation derived earlier by Lewenstein et al. [Phys. Rev. A 
51 (1995) 4617] by means of decoherence and ergodicity arguments. We show 
analytically that the thermal equilibrium state is a stationary solution of our 
rate equation. We present analytical results for the rate coefficients which are 
needed to solve the rate equations, and we give approximate formulae that 
permit their computation in practice. We solve the two sets of rate equations 
numerically and compare the results. The cooling times obtained in both 
approaches agree very well. The equilibration rates show fair agreement. 



I. INTRODUCTION 

The cooling of atoms in traps is an important tool in the study of the behavior of systems 
of Bosons and Fermions at low temperatures. Usually, the last step in the cooling process is 
evaporative cooling. For this step to be efficient, thermodynamic equilibrium of the cooled 
gas must be (nearly) attained at all times. This condition is met if the atoms interact 
sufficiently strongly. There is a number of gases in which the interaction is too weak for 
evaporative cooling to work. In such cases, one resorts to "sympathetic cooling": Another 
gas is cooled evaporatively. This gas acts as the cooling agent for the gas which is to be 
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cooled. It is usually legitimate to assume that the number of atoms in the cooling agent is 
very large, or that another mechanism is operative which keeps the cooling agent at fixed 
temperature. In either case, the cooling agent acts as a heat bath of fixed temperature. 

The method of sympathetic cooling was proposed almost twenty years ago and has 
since found widespread application (see Refs. f^J^). Exciting recent applications include 
the sympathetic cooling of ^Li Fermions in a bath of ^Li Bosons P,^], and the production 
of dual Bose-Einstein condensates with sympathetic cooling However, to the best of 

our knowledge, there does not exist until now a practicable theoretical description of this 
process. This gap is remarkable because theory would be expected to make predictions 
on the dependence of the cooling rate on various parameters defining both the interaction 
between atoms in the cooling agent and in the cooled gas, and on the trap potential, and 
might thus be helpful in improving the cooling process. We are aware of only two papers 
which deal with sympathetic cooling theoretically. Lewenstein et al. derived a master 
equation for sympathetic cooling. Unfortunately, that master equation is too complex to 
be useful, see below. Geist et al. |jT^ have formulated sympathetic cooling of Fermi gases 
in terms of a quantum Boltzmann equation. The tremendous simplification of the collision 
matrix elements achieved in this way results, however, mainly in a qualitative description of 
the cooling process. 

It is the purpose of the present work to derive rate equations for the cooling process which 
can be implemented and used practically. We do so by simplifying the master equation of 
Ref. and by reducing it to a set of rate equations. We demonstrate the utility of our 
approach by presenting analytical and numerical results which apply to the cooling process. 
In the present Section, we begin by listing the assumptions and collecting the results of 
Ref. 0. We use the notation of Ref. 0. 

The cooling gas (referred to as system B) consists of Bosons in thermal equilibrium at 
temperature Tb- We use Pb = {kTb)^^ where k is Boltzmann's constant. It is assumed 
that T does not change during the cooling process. It is also assumed that the atoms in 
system B are much heavier than those of the cooled gas. The single-particle states of 
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system B are then approximated by plane waves with energies e{k) = {hkY/{2M) with k 
the momentum and k = |A;|, with corresponding creation and annihilation operators h\k) 
and &(fc), respectively. Summations over k are always replaced by integrations. The number 
n{k) of particles with wave number k is given by 

^ '~ l-zeM-l3Be{k)] ^ ^ 

where z denotes the fugacity. While we retain Eq. (|I]) in our analytical work, in the numerical 
work of Section [V 11| below Tq is taken to be so large that the coefficients n{k) can be replaced 



by a Boltzmann distribution, 

n{k) = nBA|exp[-/5B£(A;)] , (2) 

with 

UB = (271)-^ J d^k n{k) (3) 

the density of atoms in B and = {27111^ (3 b/ MY /"^ the thermal de-Broglie wave length 
of the atoms in system B. Again, Eq. (|^) applies if the mass of the atoms in system B is 
sufficiently large. 

The atoms subject to sympathetic cooling (referred to as system A) are confined by a 
harmonic-oscillator potential with trap frequency v. We simplify the notation by confining 
ourselves to isotropic traps in three dimensions. The generalization to other dimensions, 
and to non-isotropic traps, is straightforward. The single-particle states in the harmonic- 
oscillator potential have quantum numbers n = rix, riy, n^. The integers n^, riy, are positive 
or zero and count the number of harmonic-oscillator quanta in the x, y and z directions of 
a Cartesian coordinate system. The single-particle eigenfunctions are labelled ipfi{x), the 
eigenvalues e^. The corresponding creation and annihilation operators are denoted by at 
and a^, respectively. Unless stated otherwise, we consider the atoms in system A to be 
Bosons. 

The interaction H^^b between atoms in system A and those in system B is described 
by the Hamiltonian 



Ha-b = E / dkdk'jnMk, k')alan>b^k)b{k') . (4) 

n,n' 

The matrix elements of the interaction are denoted by 'yH,n'{k, k'). These matrix elements 
contain the two-body interaction which is approximated by a delta function of strength 

C = ^ (5) 

where a is the scattering length and fi the reduced mass. Then, the matrix elements 
jri,n'{k, k') are given by 

-fH,n'{k, k') = — — / d^a; exp -i (k - k'^ x . (6) 



(27r 

The rate coefficients T^^'^ appearing in the master equation below are given by 

= dr / dWk'^^,^,{k,k')^rn,rnik',k) 

xn{k) [n{k') + 1] exp [i {e{k) - e{k') + afiu) r/Ti] . (7) 

The integer a is defined by 

" = ("^x + + T^'z) - i^x + my + m^) = {n,j. + + n^) - (n^ + n'y + n'J . (8) 

Obviously, a can be positive, negative, or zero. The last relation in Eq. (|]) restricts the 
possible values of m, m', n, n'. 

The master equation describes the dependence on time t of the reduced density matrix 
PA{t) for system A, obtained by tracing the total density matrix for systems A plus B over 
system B. Under the assumptions that the process is Markovian and that the correlation 
time for the interaction between systems A and B is much shorter than the cooling time, 
and with the help of a rotating-wave approximation, the master equation takes the form 



dt h 



Ha + -f^A-A) P{t)A 



+ CpA . (9) 



Here, Ha is the sum of the single-particle Hamiltonians for atoms in system A containing 
the kinetic energy and the harmonic-oscillator potential of the trap, while i^^-A contains 



the interaction between the atoms in system A. The action of the Liouvillean C on the 
reduced density matrix Pa(^) is given by 

n,n' ,rn,m' 

-pA{t)aiafi'alarA'^ . (10) 

Eqs. (|9|, [T0|) are not useful as they stand. Even if the accessible single-particle states of the 
harmonic oscillator are restricted to the lowest SOhu or so, the dimension of the resulting 
matrix representation of pa is enormous. Moreover, the numerical calculation of each of 
the rate coefficients T^^T presents formidable difficulties. These are compounded by the 
sheer number of such coefficients needed in the calculation. Thus, it is necessary to simplify 
Eqs. (p|,p!OD. We do so in a sequence of steps. In Section 0, we show that the decoherence 
time Tdec for the system A is inversely proportional to y^rTs where is proportional to the 
number of particles in system B, see Eq. (j^). Therefore, Tdec is the smallest time scale in the 
problem. As a consequence, the off-diagonal elements of pA are damped out very quickly, 
and Pa becomes diagonal in energy representation. The ensuing reduction of the number 
of equations and of rate coefficients is not yet sufficient, however, to lead to a practicable 
problem. We aim at a sufficiently simple rate equation. In Ref. it was pointed out that 
the degeneracy of the harmonic-oscillator states causes difficulties in the conversion of the 
master equation into a rate equation. We overcome this problem by assuming that the 
interaction between atoms in system A, although weak, lifts this degeneracy. We investigate 



two independent possibilities of reducing the problem further. In Section |T|, we average the 
density matrix over all many-body states having, for = 0, the same excitation energy 

(a fixed multiple of hu). The resulting rate equation connects diagonal matrix elements of 
Pa in energy representation. Alternatively, we introduce a factorization assumption for the 
density matrix pA (Section fV\ ). This procedure yields a simple and intuitively convincing 
form of the rate equation. The factorization assumptions used in Section ^ and, more 
implicitly, also in Section |ITT| are kin to a mean-field approximation. In Section ^ we check 
their validity by calculating the influence of the fluctuations on the rate equation which are 



neglected in the mean-field approach. In Section 0, we test the consistency of our scheme. 
We show that the equilibrium distribution of Bosons in the ground state of system A derived 
from our rate equation coincides with the one calculated by Scully |Tl|| from a different 
starting point. We then use our results for a numerical calculation. The rate coefficients 
are worked out in Section [V|. The rate equations are solved, and the results are discussed 



in Section [VI 1| . Particular attention is paid to a comparison between results obtained in the 



framework of Section 111 and of Section UM. Section VIII contains the conclusions. 



II. DECOHERENCE 

The interaction with the heat bath (i.e., with system B) not only cools system A but also 
induces decoherence: The off-diagonal elements (in energy representation) of the reduced 
density matrix Pa(^) decay in time. In this Section, we estimate the decoherence time T^ec 
and show that r^ec is proportional to l/y^rTe and, thus, much shorter than all other time 
scales in the problem. 



The common procedure |[T^ consists in calculating the time dependence of the linear 
entropy 5yi(t) defined by 

6A{t) = l-tTA{[pA{t)]') (11) 

in terms of a power-series expansion in t, assuming that at time t = the total density matrix 
p(0) of systems A + B obeys p(0) = pa(0) ® Pb(0) with tr^pA(O) = tr^(pA(0))^ = 1 and 
trBPB(O) = 1. The decoherence time T^ec is given by the inverse of the coefficient multiplying 
t or, should this term vanish, by the square root of the inverse of the coefficient multiplying 
t^. We face the second alternative because in Ref. it is assumed that tisipBHA-B) = 0. 
This causes the term linear in t to vanish. 

Second-order perturbation theory with respect to Ha-b yields 

5A{t) = -2 / dt' dt" 

n Jo Jo 

xtrA(pA(0)trB [[HA-B{t'), [^A-B(t"), Pa(0) ® Pb(0)] ] )). (12) 



The tilde indicates that the operator is taken in the interaction representation. We focus 
attention on the trace over B and consider 



X = tiB {[HA-B{t'), [HA-B{t"),PA{0)^pB{0) 



(13) 



With Hq = Ha + Hb the sum of the Hamiltonians for the atoms in systems A and B, X is 
exphcitly given by 



X=J2J2 J d'fci / d% J d^fca / d%^H,n[{ki,k[)^n,n',{hk'2) Y (14) 



where 



Y = tiBy exp{iHot' /h)ai^afi[bibp^ exp{-iHot' /h), 



We define 



•A.Hin[(t) = exp{iHAt/h)ai^afi'^ exp{-iHAt/h) 
and correspondingly for Bf:_^p{t). A straightforward calculation yields 



Y 



Ari,ri',{t'),Ari,n',{t")pA{0)\ tlB (f) (^'OPB (0) 

A^n;(t'),PA(0)A^,„- (t")] tr^ {B^,j,it')pBiO)Bj;^^,,it" 



The two traces over B yield 



(15) 



(16) 



(17) 



trB(%,^, (t')%,l(t")PB(0)) = n{h)[n{k,) + 1] 

xexp{t{eiki)-eik2))it' -t") /h) , 

trB(%,^, (t')PB(0)%£, (t")) = n{k2)[n{k^) + 1] 

X exp(i {e{ki) - 6{k2)) {f - t") /h) . (18) 

Inserting this result back into Eq. (|T^ and the latter into Eq. (^^, performing the trace 
over A, carrying out the time integrations, and expanding in powers of t, we obtain to lowest 
non-vanishing order 



6A{t) = { — ) (19) 



where 



— ) = i II f d^h f d^k2\jr{^n2{kl,k2)\' 
"dec / n J J 



X (^n{ki)[n{k2) + l]n(ni) [72(^2) + 1] 

+n{k2) [n{ki) + l]n{n2) [n{ni) + 1]) . (20) 

Here, n(n) = tr^(ata^pA(0)). 

The terms containing the occupation numbers n{k) in Eq. (|20| ) have the form 
n{ki)[n{k2) + 1] = n{ki)n{k2) + n{ki) and n{k2)[n{ki) + 1] = n{ki)n{k2) + n{k2). We focus 
attention on the terms hnear in the occupation numbers and observe that the /c-dependence 
of the remaining part of the integrand resides in |7nin2(^i; ^2)P- This quantity depends only 
on ki — k2, see Eq. (^. In the term hnear in n{ki) {n{k2)), we introduce the integration 
variables ki and k = ki — k2 {k2 and k = ki ~ k2, respectively). The integration over ki (k2, 
respectively) then yields ub, see Eq. (^). This shows that rdoc ~ ^1 

In order to implement decoherence, we follow Ref. and assume that the interaction 
H'a-a between atoms in system A is weak. Indeed, sympathetic cooling is used precisely 
when this condition is met and other procedures fail. To quantify this assumption, we 
observe that in the absence of a suitable basis for pa(^) is given by the many-body 

eigenstates |Mj) of the trap, i.e., of the three-dimensional harmonic oscillator. Here M 
denotes the total excitation energy M hu. The running index j labels the degenerate states. 
As we drop the condition = 0, the degeneracy of the states \Mj) is lifted, and even 

states belonging to different values M get mixed. We quantify the condition that the atoms 
in system A interact weakly by the requirement that the latter mixing is negligible. We 
diagonalize Ha + -f^_A-A label the resulting eigenstates by |M/i), the eigenvalues by 
e{M,fi). Here, M has the same meaning as before, while /i is a new running label. 

It turns out that for H'^_^ weak, the master equation has a curious feature. To see this, 
we take the matrix element of Eqs. (|9|, p!0|) between a state (Mi/ii| and a state \M2^2)- In- 
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spection shows that the master equation connects the time derivative of {Mifii\pA{t)\M2fi2) 
with the values of ((Mi + a)fi2,\pA(t)\{M2 + 0)^4) with a defined in Eq. (^. In other 
words, the master equation takes the form of sets of coupled equations, the equations in 
each set characterized by a fixed value of Mi — M2. The equations in different sets are 
not coupled with each other. In this situation, decoherence shows that the matrix elements 
{Mifii\pA{t)\M2fi2) with Ml 7^ M2 vanish rapidly, and that it suffices to consider only a sin- 
gle set, containing the diagonal block {Mpi\pA{t)\Mp2) ■ Within this block, we may again 
use decoherence to argue that of these matrix elements, only the ones with pi = p2 survive. 
This is what we assume from now on. 

The ensuing reduction in the number of equations is, of course, enormous: Let the total 
number of states \Mj) be A^. Then, the number of equations needed to determine Pa(^) 
is A^^ while the number remaining after we have used the decoherence argument is just N. 
However, this number is still too large, and a further simplification of the master equation 
is needed. 



III. MICROCANONICAL AVERAGE 

Because of decoherence, the reduced density matrix pA is diagonal in the basis of states 
\Mp). We now assume, in addition, that for fixed M, all matrix elements of pa are equal. 
This ergodicity-assumption seems physically reasonable: We expect the equilibration time 
of the entire system A to be much larger than the time it takes to equilibrate states which 
have (nearly) the same energy. Thus, 

{M,p,\pA{t)\M2P2) = hh,MAut^2 ■ (21) 

The Kronecker delta's express decoherence, and pM/d{M, Na) is the common value. We 
define d{M, Na) as the number of levels belonging to fixed M for a system having A^i 
particles. Transforming back to the basis of states |Mj), we see immediately that for fixed 
M, pA(t) is diagonal in that basis, too, with diagonal elements puif)- We introduce an 
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explicit notation for the states \Mj) in terms of the occupation numbers of the single- 
particle states of the harmonic oscillator. These are subject to the constraints = 
where Na is the number of atoms in system A, and X^m^m^^m = Mhu. The states \Mj) 
are written explicitly as |M{n^}). We observe that the pm introduced in Eq. pT| ) obey the 
sum rule EmPm = 1- 

Taking the partial trace of the master equation over the states \M{nrn}) with fixed M 
and using Eq. (|2l|) , we obtain 



^ 2^ J- n,n' ^ 



dt , , 

mm'nn' 



Pm 



jP d{M',NAj 



J2iM{n^^}\aiaH'ala^,\M{n^^}) . (22) 



rf(M,iV^)^ 

We note that both terms on the right-hand side of Eq. (|22D vanish unless we have fh = n' 
and n = fh'. Moreover, the terms with m = n cancel, so the case m = n is excluded from 
the summation. We also note that in the sum over M', all terms vanish but the one for 
which we have M' = M + a, with a defined in Eq. (^. The sum J^x I^I'^a)) (^{^xil 
equal to the projector onto the subspace of states with unperturbed energy M hu. The 
operator al^aa acting upon states with energy M'tiv can anyway only populate states in this 
subspace. Therefore, the projector can be replaced by the unit operator. As a result, the 
term multiplying pm' in Eq. (p2D becomes equal to Y,ii{M' {nfi}\a\arna}^a^\M' {ufi}) . Divided 
by d{M', Na), this term is nothing but the mean value of + 1] taken over the states 

with fixed M'. We denote this quantity by (^^[ra^ + 1]) M+a- The same argument is applied 
to the term in Eq. (|^) which multiplies pm and yields {nfi[n^ + 1])^. As a result, the 
master equation takes the form of a rate equation, 

= 2 51 ^H^rA(pM+a{nn[nrn + l])Af+a - PAf(?^n[?^m + ^) Al) ■ (23) 



It is straightforward to check that Y.M dPMit)/dt = 0. In Eq. (p3|), the rate coefficients 
r and the expectation values of the number operators are input quantities defined in the 
framework of our model, and the Pm's are the unknowns. 
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For later applications of Eq. (p3D, it is useful to realize that the expectation values 
depend only upon the single-particle energies and not on the single-particle wave functions. 
Therefore, we define 

f j,i = ^H^iTi ^jhu,e{rn) 5inv,£{n) ■ (24) 

Then, Eq. (^) takes the form 

dpM ^ min {K,K-a) _ 

~;|^=2 J2 ^j,j+a[PM+a{nj+ainj + l))M+a-PM{nj+ainj + l))M] ■ (25) 

a=—K j=max (0,—a) 

It is the advantage of this equation that the rates T can be calculated once and tabulated; 
the solution of Eq. (^) then requires a smaller number of steps than that of Eq. (pS)). 

To use Eq. (|23|), it is necessary to assign values to the microcanonical averages 
(^n[^»n + ^])m- We do so by replacing the microcanonical average by the canonical av- 
erage, determining the temperature T{E) by the thermodynamic relation 

^ ^K\nd{E,NA) (26) 



T{E) dE 

where d{E, N^) is the number of microstates of the Boson system A at energy E, a smoothed 
version of the degeneracy d{M,NA)- The microcanonical averages (n,i[nrra + l])Af and 
(^m[^n + 1])a/ are then approximated by nfi{T{E))[njf,{T{E)) + l] and nrfi(T{E))[nfi(T{E)) + 
1], respectively, with 

f3{E) = 1/{k,T{E)), and z{E) determined from the total number Na of particles in system 
A. The approximation neglects possible correlations between the occupation numbers n^j 
and rijfi in the states labelled M. It is, thus, related to a similar approximation used in 
Section ^ where such correlations will be discussed. 

To carry out the calculation, we need to determine the (energy-smoothed) number of 
accessible microstates d{E, Na) of system A. We do so for H'j^_y^ = 0. We observe that M 
may attain very large values, and that the degeneracy d{M, Na) increases rapidly (nearly 
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exponentially) with increasing M. If, for instance, the trap has 30 bound states and the 
system A consists of 10^ atoms, then M assumes all integers between zero and 3 ■ 10^. The 
very large degeneracy attained for large values of M causes the level density to increase very 
strongly with energy. This increase is expected to level off at a point where the majority of 
Bosons occupies the highest available single-particle level. 

We follow Ref. ||14| to compute d{E, Na). The number of accessible microstates is related 



to the grand canonical partition function S via 

K 



j=0 

oo 

= E z'^^Y.^-'^diE^N^)^ (28) 

Na=0 E 

where Khu is the energy of the highest single-particle orbital in the harmonic trap, E 
assumes values Mhu, and z denotes the fugacity. The degeneracy of single-particle states 



with energies Ej = jhu in the harmonic trap is gj = (j + l)(j + 2)/2. We introduce 
and rewrite 



X 



-/3hiy 



S = E ^"^^ E diM, Na) . (29) 

Na=0 M=0 

The microcanonical number of states d{E, N^) can be obtained from Eq. ( p9D by contour 
integration, 

d{E, Na) = fdx^fdz exp [-F{x, z)] , (30) 

with 

/ E \ ^ 
F{x,z) = i^— + Ijlnx + {Na + l)\nz + J2gj ln(l-zx^) . (31) 

The integrations in Eq. ( pOf ) can be performed using the saddle-point approximation. 
However, care has to be taken in the regime of Bose-Einstein-Condensation (BEC). We 
recall that BEC is reached in spherical harmonic traps at condensation temperatures 
kTc ^ OMN\^^hu, see e.g. Ref. ]T§. For A^^ = 400 Bosons one thus finds kT^ ^ QMu. 
In the presence of BEC the singularity due to the condensate must be excluded from the 
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saddle-point approximation |T^. The saddle points of Eq. (pO]) are difficult to obtain for 
high energies, i.e., energies beyond the maximum of the curve in Fig. |1| and especially for 
energies close to NaKKv. In this regime we simply invert the single-particle spectrum and 
solve the corresponding problem for low energies close to the (now highly degenerate) ground 
state. 

Fig. shows the number of microstates d{E, N a) as a function of the total energy E 
for a system of A^^ = 400 Bosons in a trap with K = 21. The function d{E,NA) grows 
exponentially for small energies and reaches a maximum; for energies beyond the maximum 
the finite total number of single-particle orbitals causes a decrease in d{E, Na)- The number 
of accessible microstates is a smooth function of energy. In the numerical computations we 
tabulate \iad{E, Na) for some energies and use these values for interpolation. 



1200 



800- 



400- 




4000 

E/nv 

FIG. 1. Entropy \nd{E, Na) as a function of total energy E for a system of Na = 400 Bosons 
in a trap with highest single-particle orbital K = 21. 

In concluding this Section, we note that all steps in the derivation leading to Eq. (^) ap- 
ply equally if system A consists of Fermions instead of Bosons. This results in the following 
changes. We replace r;,^(T(_E))[n^(T(_E)) + l] everywhere by r;,^(T(i?))[l—'n,f^(T(i?))]. More- 
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over, Eq. is replaced by n^(T(^)) = z{E) exp[-(3{E)e{m)]/[l + z{E) exp[-f3{E)e{7n)]]. 



IV. FACTORIZATION 

We turn to a different approximation to the master equation, written in the form of 
Eqs. (PJTO|). As in Section |II|, we assume that pA(t) is diagonal in the \Mfi) basis, and that 
- for fixed M - Eq. (|2T| ) applies. We recall the basis \M{n^}) introduced in Section |T|. 
We now change the notation, omit the letter M, and impose no further restrictions on the 
n^'s except for J2m '^m = Na- In this basis, the reduced density matrix is diagonal, and the 
master equation takes the form 

^({^4IPA(t)| W) = r™;^(2(W|a^anPA(t)ata™|{n4) 

-{{^k}\ (^nO'ma^rhO'npA (t) ling}) 

-{{n^}\PA{t)aia^alari\{n^})^ . (32) 

The last two terms are equal and combine to — 2n^[r;,^ + l]({r2^}|p^(t)|{n^}). Here, ra^ 
and denote the values of the occupation numbers for the single-particle states n and 
m, respectively, in the set {n^}. The first term yields +2n^[n^ + 1]({. . . (n^ — 1) . . . (n^ + 
1) . . .}|p^(t)|{. . . (n^ — 1) . . . (n^ + 1) . . .}). The notation indicates how the occupation 
numbers in the set {^k} are modified. The master equation becomes 

^(K}|pA(t)IK}) = 2Era 

X + 1]({. . . (n^ - 1) . . . (n^ + 1) . . .}\pAit)\{- . . (n^ - 1) . . . (n^ + 1) . . .}) 

-nnirini + l]{{n^}\pA{t)\{ng})^ . (33) 

We now take a partial trace of Eq. (p3D, summing over all occupation numbers n^^ of single- 
particle states m with m ^ mo, keeping the latter fixed. The corresponding partial trace of 
PA{t) is denoted by p(n^Q,t). This yields 

^p(^mo, ^) = 2 X! ^h'^ [ J2 ^rnirin + l]p(n^o, -l,nn + 1, t) 
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- nH[nrA + l]p{nrno,njf,,nH,t)\ . (34) 



The quantities p{nrfif^,nrfi, Ufi, t) are defined as partial traces over pA(t) excluding the single- 
particle states labelled mo, m, and n. It is easy to see that on the right-hand side of Eq. 
all terms cancel for which both m and n differ from mg. As a result, we find 

^P(^mo, t)=2 r?0:" (rirf,^ Yi^n + 1] p(^mo " 1, '^n + 1, t) 

- [Urho + 1] XI Vi^rno ,nH,t) 



-^rno Ylnn + 1] p{nrno,nH, t)^ ■ (35) 

Except for the assumptions stated at the beginning of this Section, all our steps have been 
exact. Unfortunately, the Eqs. (^5|) are not closed. We close them by introducing a factor- 
ization assumption for the partial traces of the reduced density matrix. We write 

P(^m, ^n, t) = p{nrh, t) p^Ufi, t) . (36) 

We discuss the implications of this assumption below and first deduce the results. For clarity, 
we replace the symbols p{n^, t) by Pffiin, t). This indicates more clearly that Pr^{n, t) is the 
density matrix for the single-particle state fh with occupation number n. Eq. ( pS]) takes the 
form 

Ap^„(n,t) = 2 Y TtxSnYin' + 1] Prn,{n - l,t) p,(n' + l,t) 
-[n+l]J2^' Prno{n,t) Pn{n,t) 

n' ) 

+2 E ([^ + 1] E ^' Vrn, {n + 1, t) p^{n' -l,t) 

-n EK + 1] Prnoin, t) Prf,{n', t) . (37) 

n' ) 

Summing Eq. (p^) over all n, we obtain d X^n Pmo (^5 ^) /^^ = 0- Because we have 
YliuPrhoiP'-,'^ = tr pA(t), this result is in keeping with the condition tr pA{t) = 1. We 
define the average occupation numbers Nrnit) in the single-particle state m at time t by 
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N,^{t) = J2np^,{n,t) . (38) 

n 

Multiplying Eq. (|37D by n and summing over n, we obtain the following closed system of 
rate equations for the iV^(t), 

^iV^„(t) = 2 Y: TtXm) (Nrnoit) + 1) 

-2 E r5;5(iv^(t) + i)iv^oW- (39) 

Eq. p9| ) constitutes the central result of this Section. It is a set of rate equations for the 
mean occupation numbers Nfi{t). The form of these equations is intuitively obvious and 
might even have been written down without much derivation. It is straightforward to check 
that Eq. ( ^9|) implies the condition dJ2n^nit)/dt = 0. Thus, the average particle number 
is conserved. 

The rate equations (|39|) have been derived under the tacit assumption that a finite 
number of levels in the harmonic-oscillator potential of the trap is available for occupation 
by atoms in system A. In actual fact, the trap is open, and atoms with energies above a 
critical energy are able to escape. This fact can most easily be incorporated into Eq. (|39D 
by restricting the values of m, n and mo to the bound states of the trap, and by adding 
on the right-hand side of Eq. ( pQ]) a loss term. This term allows for scattering into virtual 
harmonic-oscillator levels which are actually unbound and from which the atoms can escape. 
The derivation of the modified form of Eq. ( |39D is straightforward and yields 

^N^,{t) = 2 E ^trANrtim^^it) + 1) 
n^mo 

-2Er2S'A^r.oW- (40) 

The summation over rh' extends over the virtual states. 

The rate coefficients obey the important thermodynamic identity 

r^;^ = r;^^ ■ exp[-pB{ern -Sri)] . (41) 
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This follows immediately from the definition in Eq. (|^). For the time-independent equilib- 
rium solutions of Eq. (|39|), Eq. (^) implies the values 

^ l-ZAeM-^BSrn] ' ^ ^ 

with the fugacity za determined from the condition Y^m^r^ = ^a- 

It goes without saying that a parallel derivation applies if the atoms in system A are 
Fermions. The final rate equations are obtained by replacing on the right-hand side of 
Eq. (H) the terms (N^oi^) + 1) and (Ar^(t) + 1) hj {1 - Nr^,{t)) and by (1 - Ar^(t)), 
respectively. The equilibrium distribution is given by Eq. (^) with the minus sign in the 
denominator replaced by a plus sign. 

The central approximation Eq. (^) assumes that there are no correlations between the 
occupancies of the states m and n. The one correlation which must exist is total particle 
number conservation. Thus, we expect that the approximation Eq. (^6|) may fail whenever 
A^,^ or A^^ approach A^^, the total number of particles in system A. For Fermions, this can 
never happen. Thus, we focus attention on the case of Bosons. A critical situation arises 
if the number uq of Bosons in the ground state is comparable to A^^. A modification of 
our previous derivation is, therefore, necessary only for those terms in the master equa- 
tion which contain the functions p{no,n^,t). For these terms, we have previously assumed 
that I]n^ p{nQ,nfi,t) = Nf{{t) p{nQ,t) with Nf{{t) independent of uq. We now improve 
on this approximation by letting Nfi{t) depend on uq. This is accomplished by writing 
J2n p(^0; ''^n, t) = Nfi^TiQ, t) p^uq, t) and by using for Nfi{no, t) the ansatz 

with Ao(t) = Z^no ^0 p{no,t) and A''^(t) independent of uq. This approximation conserves 
particle number. It leads to the following modification of the rate equations. We define 
the variance 6Nq = J2noi^o)'^p{no,t) — (Z^no ^op(^o, ^))^ and have, omitting loss terms, for 
mo 7^ 

AAr^^(t) = 2 ^trCMmrnM + 1) 

n^mojO 
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-2rSiv,„(t) fivo(t) fi - ,,,,,, i + 1 1 . (44) 



The equation for A^o(^) has the form 

AiVo(t) = +2 ^ T%N^{t) fiVo(t) f 1 



dt ^ ,„ 



No{t){NA - 


AT / J \ \ 






A T /' / AT" 


AT" /J.^^ 






iVo(t)(iVA- 








No{t){NA- 





-2 E rS;>„(*) iv,(«) 1 - 1,,,,,, ) + 1 1 . (45) 



Comparing Eqs. ( pt^S]) with the original rate equations (^), we note the appearance of 
the correction term 6Nq/ {NQ{t){NA — No(t))). This term effectively reduces the coupling 
between the atoms in the ground state and the rest of the gas and is expected to increase the 
cooling time. Eqs. (|^,^) are not closed as they stand and require the solution of equations 
for higher moments of the occupation numbers. These, in turn, would not be closed. We 
address this problem in Eq. (^) below. First, we present a simple estimate of the correction 
term. We expect the term to be small both for A^o(^) ^ and for A^^ — NQ^t) <^ A^^. 
Thus, the correction term should attain its maximum at or near NQ{t) ^ (1/2) A^^. It is then 
reasonable to approximate the correction term by a simple smooth function of x = Nq/Na in 
the interval [0, 1], with very small values at the end points and a maximum near the middle. 
Qualitative support for this idea comes from Figures 1, 2, and 8 of Ref. [|1^. Using these 
Figures, we are led to the conclusion that the correction term does indeed approximately 
have the form just suggested, with a maximum value of the order of 1 per cent. It must be 
stressed, of course, that Ref. |TB| deals with equilibrium phenomena and not, as we do here, 



with equilibration processes. 

A quantitative evaluation of the correlation requires additional work. To be brief, we only 
sketch the derivation of the coupled equations which determine both, the mean occupation 
numbers A^m(^) and the ground-state correlators. A simple derivation consists in multiplying 
the master equation, Eqs. either with alaj or with alajataj^ and taking the trace. 



We use the assumptions introduced above. In particular, we assume that tr (alajpA) and 
tr (alajolafPA) are diagonal, and that the same assumption applies to the terms involving 
six creation and annihilation operators. Such terms result from the right-hand side of 
Eq. (plQl). For A: 7^ 0, we define the correlator 5^ by writing tr {alaoolaj^pj^it)) = No(t)Nj:(t) + 
Sj^{t). To obtain a closed set of equations, we use for ^ k ^ fh ^ the approximation 
that tr (ajaoata^at a^p^(t)) = No(t)Nj:(t)Nrfi(t) + Nj:{t)6^ + Nrfi{t)5j^. This approximation 
amounts to the neglect of all correlations not involving the ground state. As a result, we find 
that Eqs. retain their form, the terms -No{t)N^,,{t)5N^ / {No{t){NA - No{t))) and 

-No{t)Nri{t)6N^/ {No{t){NA- No{t))) being replaced by 5^,, and 5^, respectively. While 
Eqs. (|4^ , ^5D suggest that we have to determine a single function 6Nq/ (A'o(t)(A^A — ^oit))), 
the replacement just indicated shows that we must determine a set of functions 6^. The 
determining equations for 5^ with k read 



d 

d^^ 











\m^k,0 J 





-2Nj:{t) 



















\mi^kfi 



\fn.^kfl 



[iVrn(t) + 1] 



+2 (iVo(t) - N^,{t)) (r°;J - rjj) [6^, + NomM ■ (46) 

The last term in Eq. (^) is the "feeding term": For A^o(^) = 0, the homogeneous equation 
has the solution = 0. Deviations are due to non-zero occupation numbers of the ground 
state. 



V. CONSISTENCY 

In Refs. |]rT] , p!3| , the probability distribution for the ground-state occupation for a system 
of Bosons coupled to a heat bath was derived for the first time. In this Section, we show that 
our rate equations yield the same solution, although the system under consideration differs 
from that of Refs. [jrT],|T3[. This result then serves as consistency check for our derivation. 
Our assumptions are similar to those used in Refs. [[ri] , |T3| . More specifically, we assume 
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that the excited levels of system A are in thermal equilibrium with the heat bath (system 
B). This assumption is quantified below. 
We specialize Eq. (|35| ) to rho = 0, 

d 



-po(n,t) = 2 r^'^l n poin - l,t){NrMn-i - [n + 1] poin,t){NrMr 



+2 E + 1] Poi^ + l,t)[(iVnW)n+l + 1] 

-npo{n,t)[{Nn{t))n + l]Y (47) 

Here, (A^n(i))n = Z^n^ '^n p(^05 nfi, t)/po{n, t) is the expected value of A^n(t), given that there 
are n Bosons in the ground state. We define the cooling and heating coefficients 

= 2 E ^%{m))n. H^ = 2Y: K,imrn{t))n + l] • (48) 



This yields 

dpo(yi) 
dt 



-^Kn (n + l)po(^) - Kn-i npo{n - 1) 

npoin) - Hn+i (n + l)po(ri + 1)| . (49) 



Eqs. (|48|j49|) agree formally with Eq. (8) and the definitions following it of Ref. [jTT|, and 
with Eqs. (20,21,22) of Ref. |13|. 

The equilibrium solution of Eq. (|i9| ) has the form 

pM=Po{^)f{Kr-i/H, . (50) 

i=l 

Using this result for all = 0, . . . , Na and the constraint Y^n=o "^Poi"^) = '^o? we find for the 
normalization factor 

]_ Na i 

Po(0)-^ = Z = -E ^\{K,-,/H, (51) 

'^O i=o 3=1 

where Z is the partition function. 

We now assume that for n 7^ 0, the equilibrium mean occupation numbers {Nfi)n are 
given by a thermal distribution, 
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(iVn), 



(52) 



1 - Znexp[-(3B£n\ 

The fugacities determined by the condition 

Y,{NH)n = NA-n. (53) 

We use Eq. (|52D and the thermodynamic identity Eq. (^Tj) and find that Hn can be 
written in the form 



n^O 



(54) 



We assume that the temperature Tb of the heat bath is so low that 2;„ exp[— /^^e^] <^ 
1. Then, Hn = H becomes independent of n, and the normahzation condition Eq. (|53D 
imphes that z„ is approximately given by Zn = (A^a ~ n)z. The constant z depends upon 
temperature but is independent of n. Together with Eq. 
yields Kn = {Na — n)Hz sxpIPb^o]- Inserting this into Eq. 
condition (^Tp, we find 



Po{n) 



1 {z explPBSolY 
Z {NA-n)\ 



and the definition (^), this 
and using the normalization 

(55) 



with 



Na 



[zexp[PB£o]T 



\ 11- 



(56) 



The dependence of this result on n coincides with that of the last displayed equation on the 



left-hand side of page 023609 of Ref. [13|. This shows the consistency of our result. 



VI. RATE COEFFICIENTS 

We turn to the computation of the rate coefficients defined in Eq. (|^; various quantities 
used in this equation are in turn defined in Eqs. (^ to (^ and Eq. (|^). In what follows 
we assume that the bath particles can accurately be described by a Boltzmann distribution, 
see Eq. (^. This assumption simplifies the computation of the rate coefficients considerably 
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since it leads to a factorization of the integrals involved. Furthermore, it allows us to 
approximate in Eq. (|^ the factor [n{k') + 1] by unity since the occupation numbers of 
the bath states are small. We note that the r-integration in Eq. (|^ yields a 5-function 
and thus implies e{k) — e{k') + ahu = 0. We may therefore replace in Eq. e{k) by 
{e{k) + £{k') — ahv)/2. We obtain an integrand that is more symmetric in k and k' . Our 
starting point is thus 

rS = "^e^^B^nu J ^^J ^^kd'k'^nAk.k')ir^,n{k',k) 

—oo 

X exp [-pB{e{k) + £ik'))/2] exp [t {e{k) - £{k') + ahu) t/U] . (57) 
We introduce the oscillator length Iq = ( — ) and the dimensionless integration variables 



R = lok, k' = lok',r = x/lo,r' = x'/Iq, and t = tv^. We also define the dimensionless 
parameters a' = and 6 = ^Pb^^- Furthermore, we introduce dimensionless oscillator 
wave functions without changing our notation. Choosing Cartesian coordinates for positions 
r, "F and momenta k, k' leads to a factorization of the integrals, 

oo 

rJ5 = ^e^-'^ / dt exp (za't) [] Im„n,{t) ■ (58) 
The t-dependent integrals Im,n are given by 

oo 

Im,n{^) = I drdr' dndn' il)^{r)il)n{r)ipni{r')ipn{f'') 

.t 



exp 



exp 



2 

The normalization factor uj in Eq. is given by 



exp — K'){r — r')] . (59) 



327r4 ^ l?i Mm ^ ' 



Below we find that Im,n{t) is an even function in t and symmetric under exchange m ^ n. 
Thus, the rate coefficients ( |58D obviously fulfill the thermodynamic identity (|^). Most of 
the parameter dependence for the rate coefficients is simply contained in the normalization 
UJ, which obeys <^ z/ in realistic applications. The nontrivial and interesting parameters 
are 5 and a' . Below we discuss how the rate coefficients depend on these parameters. 
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The four integrations over r, r', k, k' can be done exactly. We first perform the integra- 
tions over r, r' and then the resulting Gaussian integrations over k, k' 0. We find 

_ /T -'-^il^:'") r(m + n-fc-/ + l/2) 

(1 + 5/4 + tv<5r+"-'=-'+^/^ ' ^ ^ 

with coefficients 

(—1)' \/m\ n\ . . 

Insertion of these results into Eq. (|58|) leads to an integral over t of the form 



oo 



Aa't 



Vnia')^ / dt- ■ (63) 

(l + l + f) 



This integral can be done exactly |jT6|. We find 

V;(a')=2i/^r(l/2-n) (^^\L \ Kj\a'\J5{l + 5/A)] . (64) 

Here -/^^^(x) denotes the modified Bessel function and n is a positive integer. We also need 
the value of Ki(a') for a' = 0. It is safe to take the limit a' ^ in the equation above. We 
find 

The exact integration of Vn constitutes an improvement over the approximate result given 
in Ref. We have reduced the computation of the rate coefficients to a six- fold sum 
{qj = rrij + rij — Ij — kj and pj = min {rrij, rij)) 

S Px Py Pz 



- - /7r\ f 1 / " '^"^ / 1 

r:^;;^ = 8 T^e^"^ E E E n Cm,,n„l, c^„n„k,T{q, + 1/2) ]V,+,^+,^+,^. 

(66) 



However, the numerical computation of the rate coefficients still presents a formidable prob- 
lem. This is due to the cancellations that occur in the sum in Eq. (|66D. These cancellations 
arise already in the two-fold sum in Eq. ([HT|). For min(m, n) exceeding ^ 20, the individual 
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terms of the sum vary by so many orders of magnitude that a numerical computation with 
double precision floating point variables leads to a total loss of precision in the final numerical 
result. Resorting to extended numerical precision is no practical solution since it increases 
the computation time by orders of magnitude. Such an increase is unaffordable in view of 
the fact that realistically, traps contain bound states of to the 30th harmonic-oscillator level 
or so. The number of rate coefficients needed in such a case is about [(1/6) ■ 30'^]^/2 ^ 10''. 

To circumvent this problem we go back to Eq. (^) and perform the Gaussian integrals 
over n and k! first. This yields 



oo / oo 



Im,n = '^nj'^ I dr^„(r) I / dr'ipmir') exp -j{r - r'f ipn{,r') ] tprnir) . (67) 



6 

We have used the shorthand notation 

7 = 2(5/(5' + 4^2) . (68) 

The integral in Eq. ( p7| ) can be viewed as the matrix element of the Gaussian two-body 
"interaction potential" G = exp [— 7(r — r')'] taken between a pair of two-body states, i.e., 

I^,n = 27r^ {fI\G\u) . (69) 

We have denoted the states by their quantum numbers as /i = {n,m) and z7 = {m,n). We 
recall that we are particularly interested in the values of these matrix elements for large 
values of the integers m, n ^ 1. In the limit of large quantum numbers a semiclassical 



evaluation of the matrix element is promising. Following Ref. [0] we use the semiclassical 
approximation and obtain 



27r 2n 

(^|G|z7) ^j^Jd^J M' e-^(^-''y^ G ( ) . (70) 



Here, = u and J^^ = fl denote the vectors (J, J') of classical actions J of the corresponding 
initial and final quantum states, respectively, and i) = (i?, i?'). We use the harmonic- 
oscillator relations r = v^2Jcos'(9, r' = \/2J' cosd' between positions and angle-action 
variables and obtain 
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2-K 2-K 

^ 

Transformation to coordinates = 0' = {d+d') /2 allows us to perform one integration 

analytically ||T^. We arrive at 

^ ' 

Here /q denotes the modified Bessel function of zeroth order. For the computation of the rate 
coefficients we perform the remaining integration over and, subsequently, the integration 
over t, numerically. 

Using Mathematica we checked that for m,n ^ 1, Eqs. ( |69D and (0) are in excellent 
agreement with Eq. (|6lD. This is exactly the regime we are interested in. As expected, the 
semiclassical approximation becomes inaccurate for \m — n\ ^ m + n. This is of no concern 
to us, however, since our exact expressions (^) and (|66D are easily computed numerically 
in this regime. We conclude that our results permit us to compute the rate coefficients 
efficiently and accurately for rather large systems. 

Fig. shows the rate coefficients for ^^Na atoms in a ^^Rb bath at temperature I/Pb = 
Thu as a function of the energy transfer (e(m) — e{n))/%v. The plotted line is the average 
value and the error bars indicate maximal and minimal rate coefficients at given energy 
transfer. In this example we choose the highest single-particle orbit at excitation energy Khu 
with = 21. The rate coefficients attain maximum values close to zero energy transfer. This 
fact supports the approximation introduced in Eq. (0). We note that the distribution of the 
rate coefficients at fixed energy transfer displays a large variance. The number of individual 
rate coefficients increases dramatically with decreasing modulus of energy transfer. The 
asymmetry between transitions with negative energy transfer (cooling) and positive energy 
transfer (heating) is due to the thermodynamic identity (|^). Approximately 2% of the rate 
coefficients had to be computed using the semiclassical method. This fraction increases with 
increasing trap cutoff K. 
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FIG. 2. Distribution of scaled rate coefficients r™^^/w for ^'^Na atoms in a ^''Rb batli at temper- 
ature = Ifiv as a function of energy transfer (e(m) — e{n)) /fiv. Average value (line), maximal 
and minimal values (error bars). 



For applications of our results, the dependence of the coohng rate on the system pa- 
rameters is of central interest. The main dependence of the rate coefficients on the system 
parameters is given by the common factor uj defined in Eq. (^). This quantity sets the 
overall time scale. As expected, uj is linear in both the density of bath atoms and the 
oscillator frequency and quadratic in the scattering length. It is a symmetric function of 
the masses of the atoms in systems A and S, attaining maximum values when one of the 
two masses is much larger than the other. The dependende of the rate coefficients on the 
parameters 5 = {m/ M)(3Bhv and a' = {M/m)a is implicit. Eq. (^) expresses the rate 
coefficients as Fourier transforms with respect to t of a three-fold product of functions Im,n- 
These functions, given in Eq. (^) and (for the semiclassical approximation) in Eq. (^) and 
Eq. ([72|), depend on t and on the parameter S. We note that Im,n is an even function of t; we 
restrict the discussion to non-negative arguments. We found that the Im,n^ are positive and 
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monotonically decreasing functions of t which vanish asymptotically for t oo. For fixed 
5 and fixed sum m + n the functions Jm,n increase with decreasing \m — n\. This behavior 
is unaffected by the integration over t and is clearly reflected in Fig. ^ For fixed 5 and 
fixed energy transfer |m — r;,| but increasing values of m + n, the values of the functions Jm,n 
increase at the origin but fall off more quickly with increasing t. In realistic applications 
we have (SbTt^v <^ 1, and this results in 5 ^ 1 unless m ^ M. For fixed m, n, the function 
Im,n increases with decreasing 5. However, this increase does not translate directly into 
an increase of the rate coefficients. For large energy transfer \m — n\ the increase is over- 
compensated by the highly oscillating exponential in the integrating over t, and the cooling 
process evolves mainly over transitions to levels that are close in energy. We also observe 
that an increase of the ration M/ m decreases 5 but increases a' and, thus, the frequency of 
the oscillations of the exponential. 

We computed the distribution of rate coefficients also for a system of ^''Rb atoms in a 
bath of ^^Rb atoms at temperature 1/ (3b = Thu. (One may assume that system and bath 
atoms are in different hyperfine states.) When comparing the distribution to the case of 
^^Na atoms in a ^^Rb bath depicted in Fig. ^ we find the following: At small energy transfer 
the Rb-Rb distribution has average and maximal values that are about three times smaller 
than for Na-Rb. However, the averages and maximal values of the Rb-Rb system decrease 
less fast with increasing modulus of the energy transfer and are similar to those of Na-Rb 
at maximal energy transfer. 

VII. NUMERICAL RESULTS 

In this Section we study the cooling process for Bosons by solving the rate equations (p5D 
and (^) numerically. We compare the results. As an example we take system A to be 
composed of Na = 400 ^'^Na atoms, while the ^''Rb bath has the temperature 1/ (3b = 7hi>. 
This temperature is slightly above the condensation temperature kTc = 6.93/iz/ for harmonic 
traps. We assume that the highest bound state in the harmonic trap has energy Khu where 
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K = 21. 

In the micro canonical approach, the rate equations are given by Eq. (^). The compu- 
tation of the input parameters {rijirij^a + 1))*/ and fj_j is described in Section |TT| and 
Section [V^, respectively. We checked that the computed expectation values {nj)M ful- 
fill the sum rules |(j + + 2)(n,)M = A^^ and + + 2)(n,)M = M 
to 1% accuracy. The homogeneous system of linear equations given by Eq. (^51) can 
be put into matrix form. We introduce the sparse matrix A which has 2K + 1 non- 
vanishing elements in each row and column and write Eq. (^) somewhat symbolically as 
dpAf/dt = J2n ^m,nPn- We note, however, that in our example the matrix A has dimension 
8401. With max (M - K,0) < N < min (M + K, KNa), the off-diagonal elements of A are 
given by 

min{K,K+M~N) 

Am,n = 2 ^ Tjj^i^_M {nj^N_M{nj + , (73) 

j=max (0,Af-Af) 

and the diagonal elements are 

min {K,M) min {K,K -a) 

Am,m = -2 ^j,j+a{nj+a{nj + 1))m • (74) 

Q=ma,x{—K,M—KNyi) j=max(0,— a) 

The limits on the summations result either from Eq. (^) or are a consequence of the iden- 
tities {riiirij + 1))a/ = for i > M and {riiirij + 1)) k{Na-i)+i = for i < /. The conservation 
of probability in Eq. (p5D is manifest since the elements in each column of A sum up to 
zero. The complete diagonalization of the matrix A is expensive and unnecessary. Instead, 
we compute only the eigenvalues with the largest real parts. For stability reasons, there 
are no eigenvalues with positive real parts. The equilibrium solution is determined by the 
zero eigenvalue. The equilibration rate is thus equal to the modulus of the real part of the 
eigenvalue with largest negative real part; the corresponding mode is damped out last. The 



eigenvalues of interest are computed using the sparse matrix solver ARPACK |T8[. Within 
our numerical accuracy we found one zero eigenvalue and no eigenvalues with positive real 
parts. For the equilibration rate we found 7eq ~ 2.7 x lO^cj, with uo defined in Eq. (0). The 
equilibrium energy E'eq = 3961.3 can easily be computed from the eigenvector belonging to 
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eigenvalue zero; alternatively, it can be obtained from the probability distribution for large 
times in a numerical integration of the rate equations (|25|). 

It is also interesting to follow the time evolution of the total energy E{t) = 
^i^J2m MpM{t) for a system of Na atoms in a Rb bath. As the initial condition we take 
Pkna — 1 vanishing values for all other probabilities. Fig. ^ shows the time evolution 
of E{t) as obtained upon integration of the rate equations. Initially, the decay is fast but 
non-exponential. This is due to the fact that many eigenvalues of the matrix A contribute 
to the cooling process. At later times the equilibrium value Eeq is approached exponentially 
fast with the rate 7eq given above. Inspection of the data used in Fig. |l] shows that the 
temperature corresponding to the equilibrium energy E^q is kT = T.Ohu. This agrees well 
with the temperature of the bath. 
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FIG. 3. Energy difference E{t) — Eeq as a function of time within the microcanonical approach 
in the Na-Rb system. 

We turn to the description of the cooling process in terms of the rate equations (p9D. 
The number of coupled nonlinear differential equations defined by Eq. (|39D is given by the 
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number of single-particle states in the trap and equals 2024 in our example. Again, we 
consider Na atoms in a Rb bath. As the initial condition we take the A^^ = 400 Bosons to 
be equally distributed over the degenerate single-particle orbitals with energy Khv while all 
other orbitals are empty. This initial condition corresponds to the situation discussed above 
for the microcanonical approach. Fig. ^ shows a plot of the total energy as a function of time. 
Initially, the decay is fast but nonexponential and can barely be distinguished from the decay 
within the microcanonical approach (compare with Fig. |^). At later times the equilibrium 
-Eeq ~ 3901 is approached exponentially fast with the equilibration rate 7eq = 1.6 x lO'^cj. 
This rate was determined from the time evolution shown in Fig. |. Alternatively, the rate 
might be obtained upon linearization of the rate equations (|39D around the equilibrium 
solution of Section 0. However, we did not pursue this point any further. The equilibrium 
energy is about 2% smaller than the corresponding value found in the microscopic approach. 
According to the data used in Fig. |l], we find a temperature of kT = G-dhu, slightly deviating 
from the bath temperature. We recall that the bath temperature enters the rate equations 
through the rate coefficients. Within our numerical accuracy, and especially because of the 
semiclassical approximation, we expect the relative error of the rate coefficients to be a few 
percent. In view of this fact, the agreement between the temperature of the bath and the 
one found in our calculation, is quite satisfactory. 
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FIG. 4. Energy difference E{t) — E^q as a function of time in the approach based on factorization 
in the Na-Rb system. 

Comparing the results obtained in the microcanonical approach and by using factor- 
ization, we note that the equihbration rate obtained from the microcanonical approach is 
about a factor ^1.7 larger than the corresponding rate found using factorization. The loss 
of energy at short times is, however, practically identical in both approaches. It is evident 
from Fig. ^ and Fig. ^ that most of the energy is removed from system A during the early 
period of the cooling process. Therefore, both approaches yield comparable predictions for 
the cooling time. For example, about 90% of the finally removed energy E{t = 0) — E^q 
have been transfered to the bath at a time tcooi ~ 0.6 x 10"^^^"^. 

We also note that the results obtained in the microcanonical approach are somewhat less 
sensitive to the exact values of the rate coefficients than those obtained using factorization. 
This is plausible since only sums over many rate coefficients enter the rate equations in the 
former case. 

Let us finally discuss sympathetic cooling of systems with equal masses. Note that 
the regime of similar masses lies somewhat outside the scope of this work. We recall that 
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the bath particles are assumed to be much heavier than the particles of the system. This 
difference in masses is also reflected by the use of different wave functions for the particles 
of the bath and the particles of the system. A detailed computation with identical wave 
functions for system and bath may thus yield quantitatively different results. With these 
cautionary remarks in mind we turn to sympathetic cooling of ^^Rb atoms in a bath of ^''Rb 
atoms. (System and bath atoms are assumed to be in different hyperfine states). The results 
apply, of course, equally to all other cases where cooled atoms and bath particles have equal 
masses. Within the microcanonical approach we found equilibration to a final state with 
energy E'eq = 3965. This corresponds to the temperature kT = T.Ohu and agrees well 
with the temperature of the bath. The equilibration rate is 7eq = 1.1 x lO^tu and the cooling 
time tcooi ~ 1-5 X lO^^u^^. Using the approach based on factorization we found for the energy 
of the equilibrated system E^q ^ SdOOhu corresponding to a temperature kT ^ G.dhu. The 
equilibration rate was 7cq ~ 0.6 x W^u; the cooling time tcooi ~ 1-5 x 10~^ct;~^. At the initial 
stages of the cooling process there is again barely any difference between the microcanonical 
approach and the one based on factorization. Nevertheless, equilibration rates vary by a 
factor of about 1.7. Comparing the results for the Rb-Rb system with the Na-Rb system 
one thus finds that sympathetic cooling times (in units of uj~^) decrease significantly with 
decreasing mass ratio m/M between system and bath particles. In our case, the factor 
(about 2.5) is similar to the factor reported at the end of Section 



VIII. CONCLUSIONS 

We have derived rate equations for the sympathetic cooling of systems composed of 
Bosons or of Fermions. The rate equations were obtained from the master equation derived 
in Ref. in a sequence of steps. First we used the perturbatively weak interaction H'^_^ 
between the particles in the system to be cooled to lift the degeneracy of the many-body 
states in the trap. This allowed us in a second step to use the decoherence argument 
and resulted in a reduced density matrix of diagonal form. In a third step we invoked an 
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ergodicity argument and assumed that the equihbration between (quasi-) degenerate states 
is much faster than the coohng process; this assumption reduced the number of independent 
diagonal elements further. From here on we proceeded along two different routes. Within the 
microcanonical approach we obtained rate equations that govern the occupation probabilities 
of sets of (quasi-) degenerate many-body states. This results in a linear problem with sparse 
matrix. The dimension of the problem increases with particle number and trap cutoff. 
Alternatively, we traced the reduced density matrix over single-particle states. Assuming 
that occupation probabilities for different orbitals factorize we obtained a nonlinear set of 
rate equations for the mean occupation probabilities of the single-particle orbitals. The 
dimension of this problem depends only on the trap cutoff. We showed how to extend this 
approach to include correlations and particle escape from the open trap, although we did 
not yet test these extensions numerically. 

We provided several checks on our assumptions. To check the consistency of our assump- 
tions we used the rate equations and computed the probability distribution for the ground 
state occupation; the results agree with those of the literature. We further showed that 
the decoherence argument is justified. Off-diagonal elements of the reduced density matrix 
decay on a time scale that is inversely proportional to the square root of the number of bath 
particles and, thus, very short. The ergodicity argument was supported by the dependence 
of the rate coefficients on energy transfer, see Fig. ^. 

Any solution of the rate equations requires as input the rate coefficients. We derived 
analytical expressions for these coefficients. A numerical evaluation of the resulting formulas 
becomes impractical, however, for transitions between high-lying single-particle states that 
are close in energy. We solved this problem by using a semiclassical approximation. We 
note that the computation of the rate coefficients is the time-consuming part in practical 
applications of the rate equations derived in this work. We discussed the dependence of the 
rate coefficients on the parameters of the system. 

The two different rate equations obtained in this paper yield results that are in semi- 
quantitative agreement. Within both approaches, the cooling times are about the same: 
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Both descriptions of the coohng process yield almost identical short-time behavior and 
practically identical final states of the system. The equilibration rates, however, differ by a 
factor of about 1.6. It is not easy to point to the origin of this difference. We recall that 
input parameters (rate coefficients) were computed to a relative accuracy of about a few 
percent. We can only speculate that the long-time behavior of the solutions to the rate 
equations are sensitive to such details. 

In practice, the choice between both approaches depends on the problem under consider- 
ation. One has to solve either a large linear problem whose dimension depends on trap cutoff 
and particle number, or a smaller nonlinear one with a dimension that only depends on the 
trap cutoff. The treatment of large systems requires the computation of a large number of 
rate coefficients. Sums of rate coefficients enter the rate equations in the microcanonical 
approach. Larger systems may become accessible more easily once these sums can be ob- 
tained without computing all terms individually. Work along this direction is in progress. 
Likewise, we have not tested yet the validity of the factorization assumption by solving the 
combined system of equations, Eqs. (^) and (^,^Sp. 

We have assumed throughout that the interaction between particles in the cooled 

system is weak. It is known |]T3[, however, that no matter how small H'^_^ is, this interaction 



cannot be neglected in the condensed state once the number of condensed atoms is sufficiently 
large. Does this statement - and the corresponding consideration for the superconducting 
state of Fermions - seriously limite our approach? We believe not. More precisely: There 
exists a generalization of our approach which overcomes the problem. It is based on the 
observation that it is sufficient to deal with the system in mean-field approximation since 
thermodynamic properties are unaffected by collective excitations . Our rate coefficients 



are defined in terms of single-particle states and single-particle energies. We may take 
these as solutions of mean-field equations. Thus, we can extend our scheme as follows. 
We follow the cooling process as described by our rate equations to the point where the 
occupation number tiq of the ground state exceeds the critical value at which interactions 
become important. We define a new set of single-particle wave functions and energies self- 
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consistently for that value of no and calculate the rate coefficients. In this way we may 
proceed to arbitrarily large values of Uq, and similarly for a system of Fermions. Obviously, 
implementation of this scheme is practical only if we succeed in speeding up the calculation 
of the rate coefficients or of sums over these coefficients. 
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